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ABSTRACT 

In  this  paper  we  discuss  the  use  of  triangular  elements  in  the  spectral  element  method  for 
direct  simulation  of  incompressible  flow.  Triangles  provide  much  greater  geometric  flexibility 
than  quadrilateral  elements  and  are  better  conditioned  and  more  accurate  when  small  angles 
arise.  We  employ  a  family  of  tensor  product  algorithms  for  triangles,  allowing  triangular 
elements  to  be  handled  with  comparable  arithmetic  complexity  to  quadrilateral  elements. 
The  triangular  discretizations  are  applied  and  validated  on  the  Poisson  equation.  These 
discretizations  are  then  applied  to  the  incompressible  Navier-Stokes  equations  and  a  laminar 
channel  flow  solution  is  given.  These  new  triangular  spectral  elements  can  be  combined 
with  standard  quadrilateral  elements,  yielding  a  general  and  flexible  high  order  method  for 
complex  geometries  in  two  dimensions.  The  natural  generalization  to  tetrahedral  elements 
in  three  dimensions  will  be  described  in  a  future  work. 
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1.  INTRODUCTION 


The  spectral  element  method  is  a  relatively  new  method  which  provides  both  the  geomet¬ 
ric  flexibility  of  finite  element  methods  and  the  exponential  accuracy  of  spectral  methods1,2,3. 
It  is  currently  being  used  for  direct  simulation  of  incompressible  fluid  flow.  With  this  method, 
the  full  Navier-Stokes  equations  can  be  solved  in  relatively  complex  geometries.  Our  goal  is 
to  increase  the  flexibility  of  the  method,  to  enable  high  accuracy  direct  simulation  of  flows 
in  complex  geometries  to  be  done  efficiently. 

The  basic  idea  of  the  spectral  element  method  is  to  use  a  tensor  product  basis  of  orthog¬ 
onal  polynomials  on  fairly  large  curvilinear  quadrilateral  subdomains  mapped  to  the  square 
[— 1,  l]2.  The  biggest  difference  between  this  approach  and  standard  finite  elements  is  that 
with  spectral  elements,  the  elements  are  much  larger  and  the  number  of  elements  is  smaller, 
allowing  the  use  of  high  order  basis  functions  (polynomials  of  order  five  to  fifteen)  on  each 
element.  Further,  the  stiffness  matrix  quadratures  are  done  “on  the  fly”  during  residual 
calculations  as  part  of  the  solution  procedure.  This  can  be  quite  efficient,  through  the  use 
of  the  tensor  product  representation  for  polynomials. 

The  spectral  element  method  provides  the  same  exponential  convergence  to  the  solution 
as  spectral  or  pseudo-spectral  discretizations.  However,  since  the  grid  in  the  spectral  element 
method  can  be  formed  by  joining  together  any  number  of  mapped  quadrilateral  subdomains, 
this  approach  provides  far  greater  geometric  flexibility  than  use  of  a  single  spectral  domain. 
With  the  spectral  method,  only  smooth  mappings  of  the  entire  domain  are  admissible  without 
sacrificing  the  convergence  rate.  The  spectral  element  method  represents  a  liberation  from 
this  restriction,  providing  great  geometric  flexibility  without  sacrificing  spectral  convergence 
to  the  true  solution. 

However,  while  the  spectral  element  method  provides  much  greater  geometric  flexibility 
than  spectral  methods,  the  restriction  to  quadrilateral  spectral  elements  does  limit  geometric 
flexibility.  Quadrilateral  domains  are  awkward,  especially  in  regions  of  high  curvature,  which 
are  typically  also  regions  requiring  fine  resolution.  Further,  irregular  geometries  such  as  rough 
walls  of  irregular  shape  are  not  easily  gridded  by  quadrilateral  elements.  Triangulation,  used 
in  unstructured  finite  volume  and  finite  element  methods,  (see,  for  example  4),  represents  a 
far  easier  method  to  create  grids  in  such  cases. 

In  addition  to  convenience  and  geometric  flexibility,  triangular  elements  have  a  funda¬ 
mental  advantage  over  quadrilateral  elements.  It  is  well  known  that  the  conditioning  of 
quadrilateral  finite  elements  degenerates  as  the  angles  at  their  vertices  approach  0  or  180 
degrees.  However  with  triangles,  in  the  limit  as  the  smallest  angle  goes  to  zero  one  retains  a 
convergent  discretization5.  This  is  not  the  case  for  quadrilaterals  using  the  standard  tensor 
product  polynomial  bases.  While  quadrilaterals  give  consistent  discretizations  if  the  min- 
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imum  vertex  angle  is  bounded  from  below,  if  the  minimum  angle  approaches  zero  as  the 
grid  is  refined,  inconsistency  occurs.  Even  when  quadrilaterals  maintain  consistency,  the 
discretization  error  using  triangular  elements  is  always  better  when  small  angles  occur. 

The  geometric  flexibility  and  favorable  approximation  properties  of  triangular  elements 
are  widely  appreciated  in  the  context  of  finite  elements.  However,  since  the  fast  tensor 
product  algorithms  used  with  quadrilateral  elements  do  not  apply  to  triangular  elements, 
triangular  spectra!  elements  do  not  seem  to  have  been  used  until  now.  A  fast  spectral 
element  method  for  triangles  and  other  regions  was  derived  previously  by  Dubiner6,  but  is 
quite  complex,  and  is  rarely  used.  Also,  Funaro7  developed  a  technique  for  treating  triangles 
for  spectral  multidomain  methods  but  does  not  advocate  its  use. 

Dubiner’s  approach  is  based  on  construction  of  an  orthonormal  basis  for  the  polynomials 
on  a  triangle.  He  shows  that  there  are  orthonormal  bases  that  are  both  well  conditioned, 
and  have  sparse  stiffness  matrices,  leading  to  fast  algorithms.  However,  to  provide  an  easy 
means  of  enforcing  inter-element  continuity,  his  basic  method  needs  to  be  modified.  He 
thus  constructs  a  new  basis  consisting  of  interior  shape  functions,  vanishing  on  the  element 
boundary,  together  with  boundary  functions  used  to  enforce  C°  interelement  continuity. 

Our  approach  and  Dubiner’s  both  yield  compution  cost  0(N3)  to  evaluate  the  residual 
on  a  triangle  with  ( N  +  l)N/2  degrees  of  freedom,  and  both  appear  relatively  stable  and 
well  conditioned.  The  central  advantage  of  our  approach  is  its  comparative  simplicity.  In 
particular,  incorporation  of  the  new  triangular  elements  into  existing  spectral  element  codes 
is  relatively  straightforward. 

The  remainder  of  the  paper  is  organized  as  follows.  In  the  next  section  we  describe  the 
design  of  the  triangular  elements  to  be  used  in  conjunction  with  our  spectral  discretization. 
Section  3  outlines  the  residual  calculation  with  fast  tensor  products  on  the  triangular  ele¬ 
ments.  An  example  calculation  is  then  given  in  section  4  to  show  that  spectral  accuracy  is 
retained  on  the  new  triangular  elements.  Section  5  describes  the  implementation  of  trian¬ 
gular  elements  in  the  solution  of  the  incompressible  Navier-Stokes  calculations,  and  a  fluid 
flow  example  follows. 

2.  DESIGN  OF  TRIANGULAR  ELEMENTS 

The  spectral  element  method  is  a  high  order  generalization  of  finite  elements.  The  stan¬ 
dard  C°  quadrilateral  finite  element  method  takes  as  trial  space  the  tensor  product  space 
QN  of  polynomials  with  terms  of  degree  at  most  N  in  either  x  or  y.  For  triangular  elements, 
one  takes  instead  the  space  PN  of  polynomials  in  x  and  y  of  degree  N.  Since  quadrilateral 
spectral  elements  are  based  on  the  trial  space  QN  too,  we  use  the  space  PN  for  triangular 
spectral  elements.  It  seems  to  be  the  only  natural  choice. 
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In  extending  the  spectral  element  method  to  triangular  elements,  it  is  natural  to  design 
a  nodal  element  method,  since  this  is  simpler  than  directly  enforcing  the  C°  continuity 
condition  between  elements,  and  simplifies  boundary  treatments.  Moreover,  by  choosing  a 
set  of  nodes  compatible  with  those  used  for  quadrilateral  spectral  elements,  we  can  make  a 
triangular  element  that  can  be  freely  intermixed  with  quadrilateral  elements. 

For  quadrilateral  elements,  one  defines  a  standard  element  on  the  unit  square,  as  shown 
in  Figure  la.  Curvilinear  elements  needed  in  the  grid  are  then  constructed  by  mapping  to 
this  standard  element.  The  same  approach  is  followed  with  triangles,  as  shown  in  Figure  lb: 
arbitrary  triangles  are  mapped  to  the  right  triangle  T  of  sides  [0, 1]. 

Nodal  Points 

We  take  as  nodal  points  the  tensor  product  Gauss  Lobatto  points  lying  within  this 
triangle,  as  shown  in  Figure  2.  If  the  quadrature  points  for  the  27V-order  Gauss  Lobatto 
formula  on  the  interval  [0, 1]  uses  points 

2  =  K,},"o 


we  take  as  nodes,  the  points: 

{(Ci.Ci)  I  o  <i,j  <  N,  i+j<N). 

To  denote  these  points  on  the  triangle,  we  will  use  the  notation  (C«>0  )(«,;)€  A  where 

A  =  {(i,j)  |  0  <  i,j  <  N,  i  +  j<N}. 

In  order  to  avoid  confusion  we  will  adopt  the  notation  (C»>  where 

□  =  {(*',»  |  0  <  i,j  <  N} 


(1) 

(2) 


for  the  usual  index  range. 

Let  T  be  the  unit  interval  [0, 1].  The  polynomial  basis  functions  on  T  are  chosen  to  be 
the  Lagrangian  interpolants  /i,  €  PN(T)  satisfying: 

MO) =  M 


where  Sij  is  the  Kronecker  delta.  They  have  the  explicit  representation: 


hi(z ) 


1 


(1  —  z2)L'n(z) 


N(N  +  l)LN(Ci)  z-Ci  ’ 

where  L yy  is  the  Legendre  polynomial  of  order  N  and  Q,  i  =  0, . . . ,  N  are  the  roots  of  L'N. 
The  representation  of  these  nodal  points  on  the  reference  triangle  T  is  given  in  Figure  2. 
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Compatibility  of  this  triangular  element  with  the  quadrilateral  spectral  elements  is  ap¬ 
parent  for  the  edges  lying  along  the  lines  x  —  0  and  y  =  0.  We  also  have  the  Gauss  Lobatto 
point  spacing  along  the  line  x  +  y  =  1,  so  this  edge  is  also  compatible  with  the  quadrilateral 
elements,  and  with  any  edge  of  other  triangular  elements. 

Quadrature  Points 

There  is  some  degree  of  latitude  in  the  selection  of  quadrature  points  within  the 
triangular  elements.  Any  choice  that  guarantees  accurate  integration  of  polynomials  of 
sufficient  order  will  be  adequate.  There  are  many  families  of  quadrature  formulas  designed 
specifically  for  triangles5.  However,  formulas  specific  to  triangles  do  not  have  the  regular 
pattern  of  quadrature  points  needed  for  exploitation  of  tensor  product  structure.  We  choose 
instead  to  use  the  tensor  product  Gauss  Lobatto  formulas.  While  this  approximately  doubles 
the  number  of  quadrature  points  used,  as  shown  in  Figure  3,  it  does  allow  exploitation  of 
tensor  product  structure.  In  the  case  of  high  order  elements,  exploitation  of  tensor  product 
structure  is  crucial.  Figure  3  illustrates  the  choice  of  quadrature  points  on  the  reference 
triangle  T.  They  are  the  points 

{(*?«'  1  Cj  )}(.-, j)gO 

where  the  y  coordinates  (Cj)  are  the  same  as  for  the  nodal  points  and  the  x  coordinates  (t/,-) 
are  different  for  each  j  or  constant  y  line.  On  each  j  line  the  coordinates  of  the  quadrature 
points  are 

4  =  o  -  Ci)  Ci 

which  are  just  the  usual  Gauss  Lobatto  points  on  the  square  mapped  into  our  triangular 
domain. 


3.  RESIDUAL  CALCULATION 

The  key  to  the  efficiency  of  the  spectral  element  method  is  the  existence  of  fast  tensor 
product  algorithms  for  evaluation  of  residuals.  In  this  section  we  describe  our  tensor  product 
algorithm  for  the  residual  calculation  on  triangles.  We  first  describe  the  polynomial  bases 
needed  and  outline  the  algorithm.  Then  we  describe  each  of  the  steps  involved  in  relative 
detail. 
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Polynomial  Bases  in  One  Dimension 

Let  N  be  the  order  of  the  spectral  element.  Define  the  one  dimensional  polynomial 
spaces: 

PN  =  polynomials  on  [0, 1]  of  order  N 

Pk,  k  <  N  =  polynomials  of  order  k 

The  N-th  order  Lagrangian  interpolation  polynomials  with  zeroes  at  the  Gauss  Lobatto 
points  are 

=  p€P"  I  p(Ci)  =  ««,  0  <  j  <  AT 

where  Sij  is  the  Kronecker  delta. 

We  also  need  the  Lagrange  basis  functions  for  Pk  with  nodes  at  the  first  k  of  the  Gauss 
Lobatto  points  {  ( «  }£lo 

s  p e P‘  | p(Cj)  =  o<j<k. 

Finally  define  the  polynomials 

u  =  4 

Then  we  have  the  relations: 

Qi  =  d* 

PN  =  span{  qi  }*0  =  span{  /,  },=0 

Pk  =  span{  dk  }*=0  =  span{  /,  }*L0 

For  later  use,  we  form  the  matrix  whose  elements  are  the  values  of  the  basis  functions 
{4}jtL:o  the  Gauss  Lobatto  points.  Define  the  matrix  G  having  elements: 

9ik  =  hiCi),  Vi,fc,  0  <i,k<N. 

Each  Ik  is  zero  on  {C»}?=d  >  one  at  (fc>  and  greater  than  one  at  the  remaining  Gauss  Lobatto 
points.  Thus  G  is  a  lower  triangular  matrix  with  unit  diagonal. 

Since  G  is  a  nonsingular  triangular  matrix,  its  inverse  is  easy  to  compute.  Let  H  —  {hij} 
be  the  inverse  of  G.  The  matrix  H  gives  the  transformation  from  function  values  at  the  Gauss 
Lobatto  points  to  polynomials  expressed  in  the  basis  {4}£L o-  That  is,  for  any  p  €  PN,  if 

N 

h *  p(  Cfc)> 

k= 0 
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then 


N 

P(x )  =  £a*',i(x)- 


In  particular,  since 

<?(£)  =  %  0  <j<k 

we  have 

771 

d*(*)  =  E 

J=0 


Polynomial  Bases  in  Two  Dimensions 

In  two  dimensions  there  are  two  polynomial  spaces  of  interest,  the  space  of  polynomials 
of  degree  N,  and  the  space  of  tensor  product  polynomials: 

PN(T)  =  polynomials  of  order  N  on  T 

Qn  =  tensor  product  polynomials  of  order  N 

where  T  is  the  right  triangle  of  side  one  shown  in  Figure  lb.  QN  is  the  space  of  polynomials 
of  degree  N  in  x  and  y  separately.  Thus  we  have  the  natural  embedding: 

Pn{T)  C  Qn 

The  notation  PN(T)  for  the  space  of  polynomials  in  two  dimensions  is  used  to  avoid  confusion 
with  our  notation  PN ,  for  the  space  of  polynomials  in  one  dimension. 

Two  different  sets  of  bases  functions  are  used  for  the  space  PN{T).  First,  given  our  choice 
of  the  tensor  product  Gauss  Lobatto  points  {(C»,Ci)  I  (*>.?)  6  A}  as  nodes  (see  Equation 
(1)  and  Figure  2),  we  define  “nodal”  basis  functions: 

4>a  s  P  €  Pn(T)  I  P((h(k)  =  6il6jk. 

The  corresponding  basis  for  PN(T)  is 

B*  =  {  foj  }  <«j)  6  a 

and  is  the  most  natural  basis  here,  though  it  is  awkward  for  computation.  An  alternate  set 
of  basis  functions  for  PN(T)  is  the  set  of  polynomials 

ft,  =  I,  i"-' 

with  the  corresponding  basis: 

B0  =  {  Pij  }  (»j)  e  a- 
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Finally,  for  the  tensor  product  polynomials  QN  we  define  the  basis  functions 

J  =  h  q3 

and  the  basis 

\>  =  {  V’tj  }  (i,j)  e  □ 

where  □  was  defined  in  Equation  (2). 

Outline  of  Computation 

Given  the  values  of  a  trial  function  u  at  the  nodes  of  the  triangle  T,  the  goal  is  to  evaluate 
the  stiffness  inner  products: 

{  a(u,  <f>ij)  }  e 

This  could  obviously  be  done  by  direct  quadrature  of  the  trial  solution  u  with  each  of 
these  test  functions  but  the  cost  would  be  0(N4).  There  are  (N  +  l)2  quadrature  points 
and  N2  -f  3iV  +  3  test  functions.  The  alternative  is  to  convert  from  the  basis  B $  to  a 
more  convenient  basis,  then  compute  quadratures  by  a  fast  tensor  product  algorithm.  The 
resulting  algorithm  is  more  complicated,  but  requires  only  0(N3)  operations,  as  in  the  case 
of  quadrilateral  spectral  elements. 

The  first  step  is  to  convert  the  trial  solution  u  expressed  in  the  basis  for  PN (T)  to 
the  basis  B ^  for  QN  via  the  operations 

B#  =$■  Bq  =>  B^. 

This  can  always  be  done,  since  PN  C  QN .  Next  we  evaluate  the  stiffness  inner  products: 

{  a(u,0o)  }  (iJ)  e  0. 

Finally,  these  inner  products  are  converted  back  to  the  basis  Bp  and  then  to  B 

{  a(ui1Pi,j)  =►  {  a{u,Pij)  }(»j)eA  {  a(ui  &%])  } (i,i)€A 

so  that  we  have  the  stiffness  inner  products  {  a(u,<f>ij)  }(i,j)eA  as  required. 

This  is,  in  outline,  the  residual  calculation.  The  total  residual  evaluation  may  be  sum¬ 
marized  in  graphical  format  as  in  Figure  4. 

The  total  cost  here  is  0(iV3),  since  each  step  is  a  tensor  product  operation  requiring 
0(N3)  operations.  The  constant  in  this  0(N3)  bound  is,  however,  larger  than  in  the  case  of 
quadrilateral  elements.  The  total  number  of  arithmetic  operations  is  about  double  that  of 
quadrilateral  elements  of  order  N.  Table  1  presents  timings  for  the  calculation  of  one  residual 
on  triangular  spectral  elements.  The  time  per  node  is  shown  to  be  proportional  with  the 


order 

time  to  compute 
element  residual 

nodes  in 
element 

time  per  node 
divided  by  order 

2 

WEBEM 1 

6 

.183 

4 

mESBrn^i 

15 

.139 

8 

48  ms 

45 

.133 

16 

338  ms 

153 

.139 

Table  1:  Execution  time  for  residual  calculation  on  triangular  spectral  elements 

order  of  the  polynomial.  This  verifies  that  the  calculation  requires  0(N3)  operations:  to  be 
exact  the  number  of  nodes  is  ( N  +  1)(jV  -f-  2)/2  so  the  number  of  operations  is  0(N(N  -f 
1  )(N  +  2)/2)  =  0((JV3  +  3N2  +  2N)/2). 

4.  VALIDATION 

The  triangular  spectral  element  Poisson  solver  is  validated  on  the  grid  shown  in  Figure  5. 
The  grid  contains  34  triangular  elements  of  varying  shapes.  Deformed  elements  are  mapped 
to  the  right  triangle  T  via  isoparametric  mappings.  Figure  6  represents  a  convergence  plot  for 
the  solution  of  Poisson’s  equation  V2u  =  /  on  the  grid  of  Figure  5  with  Dirichlet  boundary 
conditions  and  /  =  —  3.257T2  sin  rx  sin  |7rj/.  The  straight  line  convergence  in  the  log-linear 
plot  of  error  versus  N  the  order  of  the  polynomial  indicates  that  exponential  convergence 
has  been  achieved.  Hence  we  have  developed  suitable  triangular  elements  which  preserve  the 
exponential  accuracy  of  spectral  type  methods. 

Figure  7  illustrates  a  steady  state  heat  conduction  solution  in  a  fin  using  the  same  grid 
with  N  =  5  (Fig.  5):  V2T  =  q  where  q  =  0,  the  upper  and  side  boundaries  of  the  fin  root 
are  held  at  T\  =  1,  the  lower  fin  root  boundaries  are  linearly  varying  from  T\  =  1  to  T0  =  0, 
and  the  fin  surfaces  are  held  at  T0  =  0.  The  solution  is  shown  in  terms  of  isotherms. 


5.  NAVIER-STOKES  TRIANGULAR  SPECTRAL  ELEMENT  IMPLEMEN¬ 
TATION 

Direct  simulation  of  incompressible  flow  is  achieved  by  the  high  accuracy  solution  of  the 
full  unsteady  non-linear  Navier-Stokes  equations  using  triangular  spectral  elements.  The 
scheme  follows  standrad  spectral  element  procedures  described  in1,2,3,  hence  the  discussion 
will  remain  brief.  The  Navier-Stokes  equations 

du  _  ^  _  1  2  -> 

s  +  ,.v«  =  -v,+  -v« 
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V  •  u  =  0, 


where  u  is  the  velocity  vector,  p  is  the  pressure,  ana  R  is  the  Reynolds  number,  are  solved 
by  a  fractional  time-stepping  scheme  as  follows: 


u-  un 


At 


=  u  ■  Vu 


V2p  =  4-V  •  u 
y  At 


u  —  it 
At 


=  —  Vp 


un+l 


At 


This  scheme  permits  the  decoupling  of  the  pressure  and  velocity,  which,  conveniently,  can 
both  be  solved  by  Poisson  solvers.  The  nonlinear  convective  terms  are  treated  explicitly  by 
a  third  order  Adams- Bashforth  scheme  as  follows: 

2 

(u  •  Vu)n  =  £  aq(un~q  •  Vu"~?) 

,=o 


a0  =  23/12  o1  = -16/12  a2  =  5/12. 

And  hence  the  viscous  terms  are  solved  implicitly.  The  discretized  Poisson  equations  are 
solved  iteratively  by  a  preconditioned  conjugate  gradient  method.  The  preconditioners  are 
simply  the  diagonal  matrices. 

All  equations  are  solved  in  variational  form  since  derivatives  cannot  be  taken  on  the  nodal 
point  basis.  Hence,  for  example,  the  pressure  Poisson  solver  may  be  written  as 

V2p  =  j~t  V  •  u  =  ^  V  •  (u"  +  A t  u-  Vu) 

which  taken  term  by  term  gives: 


V2p  =>•  (y2p,<f>)  =  — (Vp,  V<^)  +  boundary  terms 


V  •  (u  •  Vu)  ==►  (V  •  (u  •  Vu),  <j>)  =  -(u  •  Vu,  <f>x)  -  (u  •  Vu,  <j>y)  +  b.t. 


where  u  —  ui  4-  vj.  Again,  all  these  calculations  are  implemented  with  isoparametric  map¬ 
pings  for  arbitrarily  shaped  triangular  elements. 
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6.  TWO-DIMENSIONAL  FLUID  FLOW  ILLUSTRATION 

The  Navier-Stokes  triangular  spectral  element  solver  is  illustrated  in  Figure  8  by  a  solu¬ 
tion  for  a  laminar  channel  flow.  The  four  element  triangular  grid  is  shown  along  with  the 
velocit”  profile.  The  velocity  profile  is  exact  (parabolic).  This  solution  was  obtained  from  an 
initial  solution  which  was  perturbed  by  10%  from  the  exact  solution  and  allowed  to  march 
in  time  until  it  stabilized  itself  to  the  exact  laminar  parabolic  velocity  solution.  While  the 
geometry  is  simple,  more  complex  geometry  solutions  will  be  available  by  the  conference. 

7.  CONCLUSION 

In  this  paper,  we  have  generalized  the  spectral  element  method  to  grids  containing  tri¬ 
angular  elements.  Triangular  elements  are  important  since  they  provide  greater  geometric 
flexibility  than  quadrilateral  elements.  They  have,  however,  not  been  used  previously  since 
efficient  tensor  product  algorithms  were  not  available  for  triangles.  The  algorithms  given 
here  remedy  this,  allowing  one  to  design  efficient  spectral  element  programs  using  triangular 
elements.  The  order  of  complexity  of  our  algorithms  for  triangular  elements  is  identical  to 
that  for  quadrilateral  elements,  though  the  constants  in  these  complexity  bounds  are  poorer. 
Thus  a  natural  approach  is  to  combine  triangular  and  quadrilateral  spectral  elements  us¬ 
ing  the  efficient  quadrilateral  elements  in  the  domain  interior,  and  switching  to  triangles  as 
needed  along  complex  boundaries  for  geometric  flexibility  in  modeling  incompressible  fluid 
flow. 
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Figure  1:  Elemental  mappings  for  a)  standard  quadrilateral  spectral  elements  and  b)  trian¬ 
gular  spectral  elements 
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Figure  2:  Nodal  points  used  on  the  triangular  spectral  elements 


Figure  3:  Quadrature  points  used  on  the  triangular  spectral  elements 
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Figure  4:  Total  residual  evaluation  procedure 
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T  ■  1.0 


Figure  7:  Isotherms  for  steady  state  heat  conduction  solution  on  fin  geometry  and  grid  of 
Fig.  5  with  N  =  5 
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Figure  8:  Channel  flow  solution  superimposed  on  the  triangular  grid  using  four  elements 
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